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Power-law distributions with various exponents are studied. We first introduce a 
simple and generic model that reproduces Zipf's law. We can regard this model both 
as the time evolution of the population of cities and that of the asset distribution. 
We show that our model is very robust against various variations. Next, we explain 
theoretically why our model reproduces Zipf's law. By considering the time-evolution 
equation of our model, we see that the essence of Zipf's law is an asymmetric random 
walk in a logarithmic scale. Finally, we extend our model by introducing an additional 
asymmetry. We show that the extended model reproduces various power-law exponents. 
By extending the theoretical argument for Zipf's law, we find a simple equation of the 
power-law exponent. 
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1. Introduction 

Zipf's law is the observation that the frequency n of the occurrence of various entities is an 
inverse power-law function n oc R , where the rank R is determined by the descending order of 
the frequency of each entity. Zipf 1 ) first made this observation for the frequency of occurrence 
of English words in literature; that is, the most frequent word "the" (rank R = 1) appears twice 
as many times as the second most frequent word "of" (rank R = 2). Since then, Zipf's law has 
been found in various fields of natural and social sciences including the population of cities 2 - 1 and 
the asset distribution of companies. 3 ' 4 ) This universality of Zipf's law, however, has not been 
well explained theoretically. To date, most of the available models have been specific to each 
problem. 2, 5,e ) One of the purposes of the present paper is to explain the mechanism of Zipf's law 
and thereby its universality. Our conclusion is that the essence of Zipf's law is a size change of 
each entity proportional to its size. (We made a brief report on this point previously. 7 )) 

In various fields of natural and social sciences, there are also many phenomena that exhibit 
power laws with other exponents, n oc R~ K with k ^ 1, including internet, 8 ) traffic flow, 9 ) eco- 
nomics, 10 ~ 12 ) fish school 13 ^ and family names. 14 ) We explain the mechanism of these power laws as 
well by extending the argument for Zipf's law. The exponent differs from unity whenever the size 
increase and the size decrease are asymmetric. 

In section 2, we introduce a simple and generic model that reproduces Zipf's law. We can regard 
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this model both as the time evolution of the population of cities and that of the asset distribution. 
In section 3, we explain theoretically why this model reproduces the power law. In section 4, we 
extend our model so as to reproduce power laws with other power-law exponents. Our explanation 
shows that the power-law of our model is very robust. 

Before introducing our model, let us make a tutorial remark on the relation among three types 
of distribution appearing in the present paper: the ranked distribution, the cumulative distribution 
function and the probability distribution function. The ranked distribution, including Zipf's law, 
is a cumulative distribution function with the horizontal and vertical axes being flipped. A ranked 
distribution 

nociT K (1) 

means that there are, say, R cities whose population is greater than or equal to n. In other words, 
the probability that the population of an arbitrarily chosen city is greater than or equal to n is 
given by 

P{> n) = §, (2) 

where N is the total number of cities. Changing the variable R to n on the right-hand side, we 
have 

P(> n) oc n- 1 ^ = n-\ (3) 

where b is the power-law exponent of the cumulative distribution function P(> n). Thus we 
obtain the relation between the power-law exponents of the cumulative distribution function and 
the ranked distribution. 

Next, the cumulative distribution function P(> x) is related to the probability distribution 
function p(x) in the form 

POO 

P(> x)= p(x)dx. (4) 

J X 

This means that 

p(x) oc x~ a & P(> x) oc x-^-V = x- b . (5) 

where a is the power-law exponent of the probability distribution function. 

Consequently we arrive at the relation among the three power-law exponents as follows: 

-=b = a-l. (6) 

K 

In the following sections, we introduce models that reproduce power laws in the form of the cu- 
mulative distribution function P(> x) oc x~ b . It corresponds to the power laws of the probability 
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distribution function and the ranked distribution in the forms p(x) oc x ~( b+1 ^ and n oc f? -1 / 6 , 
respectively. In particular, the case b = 1 is equivalent to Zipf's law in the form of the ranked 
distribution n oc Rr 1 . 

In the present paper, we claim that the origin of the power laws is diffusion in the logarithmic 
scale. For this reason, we also make a remark on the variable transformation of Eq (5). To see 
Eq. (5) in the logarithmic scale, we change the variable to be £ = logx. The probability distribution 
function p{x) is transformed to 

p(x)dx ~ x- a dx = e'^e^di = p(0<£- (7) 

In other words, we define the probability distribution function with respect to £ as 

P{0 = = e^. (8) 

It shows that the exponent of the exponential law of the probability distribution function in the 
logarithmic scale coincides with the power-law exponent of the cumulative distribution function. 

2. Model and simulation results 

In this section we introduce our model and show the results of its simulation. The model that 
we introduce here evolves as follows: First, we consider a set of positive values {x{} with N entities. 
Then we repeat the following procedures T times: 

Population-Distribution Model 

(i) Choose an entity i randomly from 1 < i < N. 

(ii) Add or subtract randomly the amount axi from the chosen entity with the probability 1/2: 

where a is a constant with < a < 1. For simplicity, we assume that the initial values of {xi} are 
all equal. This initial condition, however, is irrelevant in the long-time limit T — > oo. 

If we regard the values {x^ as the population of N cities, an explanation of this model from 
the viewpoint of the population distribution is that bigger cities tend to have greater population 
movement. During the time evolution we require the boundary condition that Xj should not be 
less than a lower bound xi.b.(> 0); When an operation such as X{ — > x« — axi < xi.b. is chosen 
in the step (ii), we cancel this operation. This boundary condition is required in order to prevent 
each entity x« from decreasing limitlessly. This is quite a plausible requirement; in realities, the 
population of cities, for example, cannot be less than xi.b. = 1. 

Figure 1 shows some of the results of the simulation of the above procedures for various values 
of a in the form of the cumulative distribution function P(> x) oc x~ b . (Here the number of the 
entities is iV = 10 4 .) We can clearly see that the distribution converges to power laws in all cases 
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Fig. 1. Results of the simulation for various values of a. The number of steps is T — 10 5 (dotted line), 10 6 (dashed 
line) and 10 (solid line). The straight line indicates the power law with the exponent 6=1, namely Zipf's law. 
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in the long-time limit T — > oo. Figure 2 shows the power-law exponents estimated from the results 
of our simulation. (We estimated the exponents by the least-squares method.) The exponent b is 




0.2 0.4 0.6 0.8 
a 



J 



Fig. 2. The power-law exponent 6 against the parameter a. The straight line shows the solution Eq. (31), while the 
squares indicate estimates from the results of our simulation. 

obviously close to unity for any values of a. We note that in Fig. 1 the parameter a is related only 
to the rapidity of the convergence to Zipf's law. Because of the relation (6), we obtain Zipf's law in 
the form of the ranked distribution n cx R^ 1 by flipping the horizontal and vertical axes of Fig. 1; 
the model thus reproduces the universality of Zipf's law as T — > oo. 

To show the robustness of the model, we modify it from the viewpoint of the asset distribution 
and still reproduce Zipf's law. Assume that a company tends to trade with a company of a similar 
size. Then we repeat the following procedures T times: 

Asset-Distribution Model 

(i) Choose an entity i randomly from 1 < i < N. 

(ii) For i > 1, we move the amount axi_\ from the (i — l)th entity to the ith entity, where a is a 
constant parameter with < a < 1. In other words, the following takes place: 



Xi-l 



(10) 



For % = 1, we increase all values X{ by the same amount x\ct/N . 



(ii) 
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(iii) Rearrange the entities in the descending order of the values Xj. 

An intuitive explanation of this modified model is that money flows from one company to a smaller 
company and the smaller company delivers products to the larger company. Figure 3 shows some 
of the results of the simulation of the above procedures for various values of a. (Here the number of 
the companies is N = 10 4 .) The cumulative distribution function again converges to P(> x) oc x -1 
in all cases. The operation for i = 1 in the step (ii) corresponds to the boundary condition Xj > xi.b. 
of the population-distribution model. 

The value of the lower bound xi.b. does not need to be fixed. Let us modify the population- 
distribution model and make xi.b. a changeable parameter. When the operation Xj — > Xi—axi < xi.b. 
is chosen in the step (ii) of the first model, we produce a random number S (0 < S < 1). For 

e -(*i.b.-*0 > s, (12) 

this operation is cancelled, but for 

e -(xL h -xi) < 5) ( 13 ) 

Xi becomes the new value of xi.b.. A broad change that exceeds the lower bound xi.b. is allowed 
with a small probability. Figure 4 shows the result of the simulation of the modified model for 
a = 0.3 and N = 10 4 . Compare it with Fig. 1; we can see the distribution in Fig. 4 behaves in the 
form of the Zipf 's law, while shifting to the negative direction. 

The value a does not need to be fixed as well. We can modify the population-distribution 
model and make a as a changeable parameter; we choose a randomly from < a < 1 every time 
between the steps (i) and (ii) in the population-distribution model. In other words, the value of a 
changes at each time step. Figure 5 shows the result of the simulation of the latter modified model 
for N = 10 4 . Compare it with Fig. 1; we can see the distribution in Fig. 5 behaves in the form of 
the Zipf's law. The results in Figs. 4 and 5 show that Zipf's law of our model holds out against a 
wide variety of modification. 

We see from the above variations that the essence of the model is the increase and the decrease 
proportional to the size of each entity; this is the key to the universality of Zipf's law. On the 
basis of the simpleness and the robustness of our model, we believe that Zipf's law in many fields 
of science has the same origin of "proportional change." 

3. Continuum Limit and Universality 

In this section we give a theoretical explanation as to why the model in the previous section 
reproduces Zipf's law. Let us consider the limit N — > oo. Hence we regard the set of variables 
{xj} as a continuous variable x. We consider in this limit the time evolution of the probability 
distribution function p(x). We note that Eq. (9) gives probability flows in the directions 

x — > x ± ax = (1 ± a)x. (14) 
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Fig. 3. Results of the simulation of the asset distribution for various values of a. The number of steps is T = 10 5 
(dotted line), 10 6 (dashed line) and 10 7 (solid line). We normalized the results by the asset of the highest entity. 
The straight line indicates the power law with the exponent 6=1, namely Zipf's law. 
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Fig. 4. The result of the simulation for a — 0.3 with a broad lower bound. The number of steps is T — 10 7 (dotted 
line), 4 x 10 7 (dashed line) and 10 8 (solid line). The straight line indicates the power law with the exponent 6 = 1, 
namely Zipf's law. 
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Fig. 5. The result of the simulation with randomly chosen a (0 < a < 1). The number of steps is T = 10 5 (dotted 
line), 10 6 (dashed line) and 10 7 (solid line). The straight line indicates the power law with the exponent 6=1, 
namely Zipf's law. 
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The time-evolution equation of the probability distribution function p(x) is therefore given by the 
following: 

^-p(x,t) = - 7 p(M) + \ ~T^~ P f T^-'O + \t~jL — ^ (Vl - '*! > ( 15 ) 
at 21 — a \ 1 — a J 21 + a \l + a / 

where p(x, t) is the probability distribution of x at time t and 7 is a constant determining the unit 

of time. The first term on the right-hand side represents a flow from the point x out to the points 

(1 ± a)x, while the second and third terms represents flows from x/(l =L a) into x. 

The coefficients 1 / (l±a) in front of p in Eq. (15) are necessary in order to satisfy the probability 

conservation. To see this, we integrate Eq. (15) over x as 

d f°° f°° 1 7 f°° fx \ 

-J o P (x,t)dx = -7^ p ( x ,t)dx + -—J o P[—,t)dx 

+ r p{t^A^- (16) 



2 1 + a J \1 + a 

Changing the variables to be 



*' = -^-, ^=-^-, (17) 
1 — a 1 + a 



we can rewrite Eq. (16) as 



— J p(x,t)dx = —7 J p(x,t)dx + — j p(x',t)dx' 

COD 

I 

+ 



lj™ P (x",t)dx" = 0. 



Thus the total probability is conserved. (In the present case, the probability would be conserved 
without the coefficients 1/(1 ± a). However, they are necessary when Eq. (15) is generalized later; 
see Eq. (35).) 

Let us see Eq. (15) in the logarithmic scale; we change the variable to be £ = log a;. The 
probability distribution function is transformed to 

p(x, t)dx = p(et, t)e?dt = p(£, i)d£. (18) 

In other words, we define 

p(S,t)= P (x,t)x = p(et,t)e*. (19) 

By rewriting p in terms of p as 

p(x,t) = -p(logx,t)=e-Zp(tt), (20) 

x 

fx \ 1 =F a _ / x \ 

P 7— >< = pi log——,* J 

\1 =F a J x \ 1 =f a J 

= (lTa)e- 5 j5(e-log(lTa),t), (21) 
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we transform the evolution equation (15) to 

§iffi,t) = -rmt) + lp(S + P-,t) + lp(Z-P+,t), (22) 

where 

(3± = ±log(l±a)>0. (23) 

We can regard the evolution equation (22) as a random walk from £ to £ ± (3±, that is, a random 
walk with asymmetrically fixed step sizes. 

Let us consider the meaning of the asymmetric random walk. We carry out the Taylor expansion 
of the second and third terms on the right-hand side of Eq. (22). Assuming O(0±) — 0, we can 
rewrite Eq. (22) as 

fiKct) = -2 7 p(e,t) +7 (m,t) + p-%m,t) + %$mtj) 

+7(p(e,t)-/? + ^(e,t) + ^^j5(e,t)). 

By rearranging the right-hand side of Eq. (24), we have 

§- t p(t, t) = 7 ^+±^|lp(£, t) + 7 (-/3+ + 0-)^m t). (25) 

The first and second terms on the right-hand side of Eq. (25) correspond to the diffusion and 
advection terms, respectively. Consequently we can regard the equation of the asymmetric random 
walk, Eq. (22), as an advection-diffusion equation (25) approximately. For < a < 1, the advection 
coefficient is 

-(3+ + 13. = - log(l + a) - log(l - a) = log — ^ > 0. (26) 

1 — Q z 

It suggests that the probability distribution function in the logarithmic scale j5(£, t) moves to the 
negative direction while diffusing. 

Let us now concentrate on finding stationary solutions p s {0 of Eq. (22): 

-2 7J 5 S (£) + 7 j5 s (£ + /?_) + 7 p s (£ - 13+) = 0. (27) 

We seek solutions of the form 

p s ex e~ bi . (28) 

(See Appendix A for the validity of this assumption.) By plugging Eq. (28) into Eq. (27), we have 

-2 7 e- b « + je- b ^e- bl3 - + 7e -^e 6/3 + = 7 e- b «(-2 + e~ bp - + e b/3 +) = 0, (29) 

or 



e bf> ++e -W- =2 . 



(30) 
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Rewriting Eq. (30) using Eq. (23), we obtain the following equation for 6: 

f(b)^ l -{l + a ) b + l -(l-a) b = l. (31) 

Figure 6 shows the function f(b) for various values of a. Though the form of f(b) changes with a, 
the two solutions of Eq. (31), 6 = and 1, are stable. Obviously b = is a solution of Eq. (31). 




Fig. 6. The function f(b) defined in Eq. (31) for various values of a. The curvature of f(b) increases as a becomes 
greater. 

The solution b = means that the stationary solution p s is a uniform distribution: 

p s oc e~ b ^ = const. (32) 

As we suggested in Eq. (25), our model can be regarded as diffusion with advection in the loga- 
rithmic scale. We hence expect that the solution b = 0, or a uniform distribution appears only in 
the case of periodic boundary conditions. For our model with the lower bound xu,., the solution 
6 = 1 survives: 

ps oc (33) 

According to the argument at the end of Sec. 1, the solution (33) means P(> x) oc x . This 
is in good agreement with our simulation data in Fig. 1. The present argument indicates that a 
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random walk with asymmetric step sizes in the logarithmic scale is the essence of the universality 
of Zipf 's law. 



4. Extension of the model 

So far, we have focused on Zipf's law (b = 1). However, various power-law exponents (b ^ 1) 
have been reported in many phenomena including economics 10-12 -* and family names. 14 ' In this 
section, we extend our model in order to reproduce power-law distributions with other exponents. 

The extension is basically modification of the asymmetry of the random walk. We introduce 
two types of extension. In both extensions, the initial and the boundary conditions are the same 
as in the model in Sec. 2. Extension A of the model evolves as follows: 

(i) Choose an entity i randomly from 1 < i < N. 

(ii) Add the amount a+Xj or subtract the amount a_Xj randomly from the chosen entity with the 
probability 1/2. 

Xi — ► Xi ± a±Xi. (34) 

Here the parameters a+ and a_ are made to be different; a + is the growth parameter and a_ is 
the anti-growth parameter. We restrict ourselves to the case < ot + < a_ < 1; otherwise, we can 
have a meaningless solution b < (See Eq. (38) below). 

Figure 7 shows some of the results of the simulation of the above procedures for various values 
of a and N = 10 5 . We can still see the power-law distribution, but the power-law exponents b 
changed from unity. The power-law exponent b takes various values in the modified model while it 
is fixed to b = 1 for any values of a in Sec. 2. 

We can easily modify the argument in Sec. 3 so as to analyze the above modified model, 
Extension A. Instead of Eq. (15) we now write the time-evolution equation of p(x) as 

d . . , . 1 7 fx \ 1 7 fx \ 

- P (x,t) = - 1P (x,t) + -— - P (j^.tj + - 2TT ^p [tt^'V ■ (35) 

Note that the probability is still conserved. By rewriting p in terms of p, we transform the evolution 
equation (35) to 

t) = - 7 p(£, t) + +p-,t)- \m -p+,t), (36) 

where 

(3± = ±log(l±a±) > 0. (37) 
Note the difference between Eq. (23) and Eq. (37). The equation for the exponent b now reads 

f(b)^\(l + a + ) b + l -(l-a^ b = l. (38) 
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Fig. 7. Results of the simulation of Extension A for various values of a. The number of steps is T — 10 6 (dotted 
line), 10 7 (dashed line) and 10 s (solid line). 
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Fig. 8. The function f(b) defined in Eq. (38) b for various values of a± with the ratio a+/a~ = 0.7 fixed. The 
curvature of /(&) increases as a± become greater. 
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Figure 8 shows the function f(b) for various values of a± with = 0.7. Though the 

solution b = is still stable, the other solution depends on a±. Owing to the asymmetry between 
a + and a_, the relevant solution of Eq. (38) now depends on a + and a_. Figure 9 shows that the 
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Fig. 9. The power-law exponent b against the parameter a_ while fixing a+/a. 
solution of Eq. (38), while the squares are the simulation results. 



0.7. The solid line shows the 



simulation results of Extension A (Fig. 7) are in good agreement with the solution of Eq. (38). 

Now we introduce Extension B of the model that also reproduces power-law exponents b ^ 1. 
Extension B of the model evolves as follows: 

(i) Choose an entity i randomly from 1 < i < N. 

(ii) Add or subtract randomly the amount axi from the chosen entity with the probabilities p+ 
and p-, respectively (p_ + p_ = 1). In other words, 



(39) 



Xi + axi with probability p + . 

with prabability p_. 

Figure 10 shows some of the results of the simulation of Extension B. We again see power laws 



with b ^ I. The time-evolution equation of p(x) now reads 
d 



dt 



p{x,t) 



-jp(x,t) + P-- p{- 



1 



a 



1 



a 



-PKTTT 'fi- 



at 



1 + a 



(40) 
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Fig. 10. Results of the simulation of Extension B for various values of a with p+ = 1/3 and p_ —2/3. The number 
of steps is T = 10 6 (dotted line), 10 7 (dashed line), 10 8 (solid line). 
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By rewriting p in terms of p, we transform the evolution equation (40) to 
d 



at 



where 

0± = ± log(l ± a) > 0. 
The equation for the exponent 6 now reads 

f(b)=p^l-a) b + p + (l + a) b = l. 



(41) 



(42) 



(43) 



Figure 11 shows the function /(ft) for various values of a with p + /p_ = 1/2. Though the 
solution 6 = is still stable, the other solution depends on a. Owing to the asymmetry between 
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Fig. 11. The function /(&) defined in Eq. (43) for various values of a with p+ = 1/3 and p_ = 2/3 fixed. The 
curvature of f(b) increase as a becomes greater. 

p+ and p-, the relevant solution of Eq. (43) depends on a, p+ and p-. Figure 12 shows that the 
simulation results of Extension B (Fig. 10) are in good agreement with the solution of Eq. (43). 



5. Summary 



In the present paper, we explained the origin of the universality of Zipf 's law theoretically and 
demonstrated it numerically. We claim that the essence of Zipf 's law is a size change of each entity 
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Fig. 12. The power-law exponent b against the parameter a while fixing p+/P- = 1/2. The solid line shows the 
solution of Eq. (43), while the squares are the simulation results. 
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proportional to its size. The simulation results of our simple and generic model suggest that the 
present explanation of Zipf's law is applicable to various phenomena both in natural and social 
sciences. 

We also extended the model slightly and reproduced power laws other than Zipf's law. We 
revealed that some asymmetry between the increase and the decrease causes the deviation from 
Zipf's law. Equations (38) and (43) may be useful in estimating parameters in actual power-law 
phenomena. 

Appendix A: Exponential Solutions 

In this appendix, we discuss the validity of the assumption (28) for the solution of Eq. (27). 
Carrying out the Fourier transform of Eq. (27), we obtain 

-2&(k) + e lk ^g{k) + e'^-fsik) = 0, (A-l) 

or 

e ik/3 + + e -ik{3_ = 2> ( A . 2 ) 

where 

Fs(k) = J mV^dt. (A3) 
By substituting Eq. (23), we can write Eq. (A- 2) in the form 

(1 + a) ik + (1 - a) ik = 2. (A-4) 

Let us look for solutions k in the complex plane by writing 

k = k r — ib, (A-5) 

where k r and b are real numbers. If we have solutions with k r = only, the assumption Eq. (28) 
is justified. Substituting Eq. (A-5) for k in Eq. (A-4), we have 

(1 + a) b e ikrl3 + + (1 - a) b e' ik ^- = 2. (A-6) 

To satisfy Eq. (A-6), we need the phase factors to be real: 

or 

_ e ikr/3+ = e -ik r P_ = L ( A .g) 

In the case of Eq. (A-7), we have 

2mir 2nir , . . 

K = — = IT > (A - 9) 
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while in the case of Eq. (A-8) we have 

P+ P- 

where m and n are integers. In order for the solutions (A-9) and (A-10) to exist, the ratio of the 
step sizes (3+ and /3_ must be a rational number. It means that the random walker has a possibility 
of returning exactly to the starting point. We consider that such cases are quite exceptional; if we 
choose a general value of a, the ratio /?+ //3_ is an irrational number in general. Hence we neglect 
the solutions (A-9) and (A-10) and assume k r = 0. Therefore the remaining solutions have the 
form k = —ib. We can thus justify the assumption (28). 
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